Overview

Purpose: Generate multi-panel volcano plots for differential expression analysis across multiple predictors (leaf tissue, phosphorus treatment, inv4m genotype), highlighting genes selected by Mahalanobis distance outlier detection.

Input: CSV file containing differential expression results with Mahalanobis outlier identification (is_selected_DEG column).

Output: - Multi-panel volcano plots (3-panel, 2-panel, and single-panel versions) - Annotated results table (CSV) - LaTeX tables of top labeled genes by predictor × regulation

Setup

Load Libraries

library(dplyr)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(ggtext)
library(kableExtra)
library(readr)

Data Import and Filtering

Purpose: Load differential expression results and filter to predictors of interest.

Approach: Import CSV and subset to leaf tissue, phosphorus treatment, and inv4m genotype effects. The input contains pre-calculated Mahalanobis distances and outlier flags from upstream analysis.

Expected outcome: Filtered dataset with three predictor levels, retaining all columns including is_selected_DEG for gene labeling.

predictor_order <- c("Leaf", "-P","Leaf:-P")

# Load DE results with Mahalanobis outlier detection
effects_df <- read.csv(file.path(paths$intermediate, "predictor_effects_leaf_interaction_model.csv")) %>%
  filter(predictor %in% predictor_order ) %>%
  droplevels()

# Set predictor order: Leaf, -P, Inv4m

effects_df$predictor <- factor(effects_df$predictor, levels = predictor_order)

# Diagnostics
cat("Predictors in dataset (ordered):", levels(effects_df$predictor), "\n")
## Predictors in dataset (ordered): Leaf -P Leaf:-P
cat("Total genes:", nrow(effects_df), "\n")
## Total genes: 72033
cat("Genes per predictor:\n")
## Genes per predictor:
print(table(effects_df$predictor))
## 
##    Leaf      -P Leaf:-P 
##   24011   24011   24011

Gene Selection for Labeling

Purpose: Identify top genes to label in volcano plots based on Mahalanobis distance ranking.

Approach: For each predictor and regulation direction, select the top N genes (default N=10) that are marked as is_selected_DEG and have locus_label annotation in selected_DEGs_curated_locus_label.csv. Load curated labels from the external file, coalesce with existing labels, then filter for genes with valid annotations to label in the plot.

Expected outcome: Data frame of genes to label, with custom label text from curated annotations.

# Load curated labels from intermediate results
curated <- read.csv(file.path(paths$intermediate, "selected_DEGs_leaf_interaction_model.csv")) %>%
  select(
    regulation,
    predictor,
    gene,
    locus_symbol,
    locus_label,
    desc_merged,
    mahalanobis,
    logFC,
    neglogP
  )

# Set predictor order for curated data
curated$predictor <- factor(curated$predictor, levels = predictor_order)

# Define number of labels per predictor and regulation
label_counts <- data.frame(
  predictor = c("Leaf", "Leaf", "-P", "-P", "Leaf:-P","Leaf:-P"),
  regulation = c("Downregulated", "Upregulated",
                 "Downregulated", "Upregulated", 
                 "Downregulated", "Upregulated"),
  n_labels = c(21, 20, 10, 15,10,10),
  stringsAsFactors = FALSE
)

# Select genes for labeling based on specified counts
to_plot <- curated %>%
  ungroup() %>%
  filter(!is.na(locus_label) & locus_label != "") %>%
  arrange(regulation, predictor, -mahalanobis) %>%
  left_join(label_counts, by = c("predictor", "regulation")) %>%
  group_by(regulation, predictor) %>%
  mutate(row_num = row_number()) %>%
  filter(row_num <= n_labels) %>%
  select(-row_num, -n_labels) %>%
  arrange(-mahalanobis) %>%
  ungroup() %>%
  mutate(predictor = factor(predictor, levels = predictor_order)) %>%  
  data.frame()

# Diagnostics
cat("Total curated genes:", nrow(curated), "\n")
## Total curated genes: 122
cat("Genes with locus_label:", sum(!is.na(curated$locus_label) & curated$locus_label != ""), "\n")
## Genes with locus_label: 73
cat("\nSelected DEGs by regulation and predictor:\n")
## 
## Selected DEGs by regulation and predictor:
with(curated, print(table(regulation, predictor)))
##                predictor
## regulation      Leaf -P Leaf:-P
##   Downregulated   18 17      13
##   Upregulated     19 17      11
cat("\nGenes selected for labeling:", nrow(to_plot), "\n")
## 
## Genes selected for labeling: 52
cat("Labels per predictor/regulation:\n")
## Labels per predictor/regulation:
with(to_plot, print(table(regulation, predictor)))
##                predictor
## regulation      Leaf -P Leaf:-P
##   Downregulated   12 10       5
##   Upregulated      6 14       5
# Show configured label counts
cat("\nConfigured label counts:\n")
## 
## Configured label counts:
print(label_counts)
##   predictor    regulation n_labels
## 1      Leaf Downregulated       21
## 2      Leaf   Upregulated       20
## 3        -P Downregulated       10
## 4        -P   Upregulated       15
## 5   Leaf:-P Downregulated       10
## 6   Leaf:-P   Upregulated       10

Visualization

Three-Panel Volcano Plot (Leaf, -P, Inv4m)

Purpose: Create three-panel volcano plot showing logFC vs. FDR across all predictors.

Approach: Use faceted ggplot with: - Points colored by regulation status (down/Bellow threshold/up) - Text labels for top Mahalanobis-selected genes - Free x-axis scales to accommodate different effect size ranges - Expression formatting for axis labels

Expected outcome: Three-panel figure saved as high-resolution PNG.

# Custom formatting function for y-axis
scaleFUN <- function(x) sprintf("%.0f", x)

# Set regulation based on high-confidence DEG flag
effects_df$regulation[!effects_df$is_hiconf_DEG] <- "Unregulated"
to_plot <- to_plot %>% filter(
  !(predictor=="Leaf" & regulation=="Upregulated" & neglogP < 2.5)
  )

effects_df$predictor_label <- effects_df$predictor
effects_df$predictor_label[effects_df$predictor=="Leaf:-P"] <- "Leaf × -P"
# Generate three-panel plot
gene_volcano_3 <- effects_df %>%
  ggplot(aes(x = logFC, y = neglogP, col = regulation)) +
  ylab(expression(-log[10](italic(FDR)))) +
  xlab(expression(log[2]("Fold Change"))) +
  scale_y_continuous(labels = scaleFUN,expand = expansion(mult = c(0, 0.1)) ) +
  geom_point(alpha = 0.3, size = 0.5) +
  geom_text_repel(
    data = to_plot,
    aes(label = locus_label),
    color = "black",
    fontface="bold.italic",
    segment.color ="grey",
    force = 2,
    size=7,
    min.segment.length = 0,
    max.overlaps = 50
  ) +
  scale_color_manual(
    name = "",
    values = c("#00AFBB", "grey", "#bb0c00"),
    labels = c("Downregulated", "Bellow threshold", "Upregulated")
  ) +
  facet_wrap(. ~ predictor, scales = "free_x" ) +
  guides(color = guide_legend(
    reverse = TRUE, 
    override.aes = list(size = 6)
  )) +
  theme_classic2(base_size = 30) +
  theme(
    plot.title = element_blank(),
    strip.text = element_blank(),
    #strip.text = element_markdown( size=35),
    legend.position = c(0.87, 0.95),
    legend.background = element_rect(fill = "transparent"),
    legend.spacing = unit(0, "line"),
    legend.box.spacing = unit(0, "line"),
    legend.text = element_text(size = 25),
    legend.direction = "vertical",
    plot.margin = margin(0, 5.5, 5.5,5.5, "pt"),
    strip.background = element_blank()
  )



# Load lipid volcano
lipid_volcano_3_file <- file.path(paths$intermediate, "lipid_volcano_3.rds")
if (file.exists(lipid_volcano_3_file)) {
  lipid_volcano_3 <- readRDS(lipid_volcano_3_file)
  cat("Loaded lipid volcano plot\n")
} else {
  warning("Lipid volcano plot not found. Run differential_lipid_analysis_leaf_treatment_interaction_model.Rmd first.")
}
## Loaded lipid volcano plot
# Load transcript MDS
transcript_MDS_file <- file.path(paths$intermediate, "transcript_MDS_12.rds")
if (file.exists(transcript_MDS_file)) {
  p_transcript_MDS <- readRDS(transcript_MDS_file)
  cat("Loaded transcript MDS plot\n")
} else {
  warning("Transcript MDS plot not found. Run differential_expression_leaf_treatment_model.Rmd first.")
}
## Loaded transcript MDS plot
lipid_MDS_file <- file.path(paths$intermediate, "lipid_MDS_23.rds")
if (file.exists(lipid_MDS_file)) {
  p_lipid_MDS <- readRDS(lipid_MDS_file)
  cat("Loaded lipid MDS plot\n")
} else {
  warning("Lipid MDS plot not found. Run differential_lipid_analysis_leaf_treatment_interaction_model.Rmd first.")
}
## Loaded lipid MDS plot
compound_MDS<- ggpubr::ggarrange(
  p_lipid_MDS,
  p_transcript_MDS,
  nrow=2,
  align ="hv")

ggsave(
  filename = file.path(paths$figures, "compound_MDS.png"),
  plot = compound_MDS,
  height = 14,
  width = 8,
  units = "in",
  dpi = 150
)

compound_volcano<- ggpubr::ggarrange(
  lipid_volcano_3,
  gene_volcano_3,
  ncol=1,
  align ="v"
  )

ggsave(
  filename = file.path(paths$figures, "compound_volcano.png"),
  plot = compound_volcano,
  height = 12,
  width = 14,
  units = "in",
  dpi = 150
)

 compound_plot<- ggpubr::ggarrange(
  compound_MDS,
  compound_volcano,
  ncol=2,
  widths =  c(0.5,1)
)

ggsave(
  filename = file.path(paths$figures, "compound_plot.png"),
  plot = compound_plot,
  height = 12,
  width = 20,
  units = "in",
  dpi = 150
)


compound_volcano

Two-Panel Volcano Plot (Leaf, -P only)

Purpose: Create two-panel volcano plot focusing on Leaf and -P predictors.

Approach: Filter data to exclude Inv4m, adjust legend position for two-panel layout.

Expected outcome: Two-panel figure saved as high-resolution PNG.

# Filter for Leaf and -P only
effects_2panel <- effects_df %>%
    mutate(predictor = factor(predictor, levels = predictor_order)) %>%
  filter(predictor %in% c("Leaf", "-P")) %>%
  droplevels()

to_plot_2panel <- to_plot%>%
    mutate(predictor = factor(predictor, levels = predictor_order)) %>%
  filter(predictor %in% c("Leaf", "-P")) %>%
  droplevels()
# to_plot_2panel%>% filter(locus_label=="hhtmp5")
# to_plot_2panel%>% filter(locus_label=="hhtmp5")
# to_plot_2panel$locus_label[to_plot_2panel$locus_label=="hhtmp5"] <- NA
# Generate two-panel plot
volcano_2panel <- effects_2panel %>%
  ggplot(aes(x = logFC, y = neglogP, col = regulation)) +
  ylab(expression(-log[10](italic(FDR)))) +
  xlab(expression(log[2]("Fold Change"))) +
  #scale_x_continuous(expand = expansion(mult = c(0.3, 0.3)),labels = scaleFUN) +
  scale_y_continuous(expand = expansion(mult = c(0., 0.2)),labels = scaleFUN) +
  geom_point(alpha = 0.3, size=0.3) +
  geom_text_repel(
    data = to_plot_2panel,
    aes(label = locus_label),
    color = "black",
    segment.color="grey",
    force = 1,
    fontface="bold.italic",
    size=5.5,
    fontface = "italic",
    min.segment.length = 0,
    max.overlaps = 50
  ) +
  scale_color_manual(
    name = "",
    values = c("#00AFBB", "grey", "#bb0c00"),
    labels = c("Downregulated", "Bellow threshold", "Upregulated")
  ) +
  facet_wrap(. ~ predictor, scales = "free_x") +
  guides(color = guide_legend(
    reverse = TRUE, 
    override.aes = list(size = 5)
  )) +
  theme_classic2(base_size = 30) +
  theme(
    # plot.title = element_blank(),
    strip.text = element_markdown(size=35),
    strip.background = element_blank(),
    axis.title.y = element_text(size=25,margin = margin(r = -10)),
    # axis.title.x = element_blank(),
    #plot.margin = margin(5.5, 5.5, 30, 5.5, "pt"),
    legend.position = c(0.15, 0.99),
    legend.background = element_rect(fill = "transparent"),
    legend.spacing = unit(0, "line"),
    legend.box.spacing = unit(0, "line"),
    legend.text = element_text(size = 20)
  )

Single-Panel Volcano Plot (Leaf:-P only)

Purpose: Create single-panel volcano plot focusing exclusively on Inv4m genotype effects.

Approach: Filter data to Inv4m only, optimize layout for single panel with centered legend.

Expected outcome: Single-panel figure saved as high-resolution PNG.

# Filter for Inv4m only
effects_1panel <- effects_df %>%
  filter(predictor == "Leaf:-P") 

to_plot_1panel <- to_plot %>%
  filter(predictor == "Leaf:-P")

# Generate single-panel plot
volcano_1panel <- effects_1panel %>%
  ggplot(aes(x = logFC, y = neglogP, col = regulation)) +
  ylab(expression(-log[10](italic(FDR)))) +
  xlab(expression(log[2]("Fold Change"))) +
  scale_y_continuous(labels = scaleFUN, breaks=0:6) +
  geom_point(alpha = 0.3, size = 2) +
  geom_text_repel(
    data = to_plot_1panel,
    aes(label = locus_label),
    segment.color ="grey",
    force = 2,
    fontface = "bold.italic",
    size=6,
    color = "black",
    min.segment.length = 0,
    max.overlaps = 50,
    size = 7
  ) +
  scale_color_manual(
    name = "",
    values = c("#00AFBB", "grey", "#bb0c00"),
    labels = c("Negative", "Bellow threshold", "Positive")
  ) +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  theme_classic2(base_size = 22) +
  theme(
    plot.title = element_blank(),
    legend.position = "top",
    legend.justification = "right",
    legend.background = element_rect(fill = "transparent"),
    legend.spacing = unit(0, "line"),
    legend.box.spacing = unit(0, "line"),
    legend.text = element_text(size = 20),
    legend.direction = "horizontal"
  )


ggsave(
  filename = file.path(paths$figures, "volcano_plot_LeafxP.png"),
  plot = volcano_1panel,
  height = 7,
  width = 7,
  units = "in",
  dpi = 150
)
cat("Figure saved: volcano_plot_LeafxP.png (7 x 7 inches @ 150 dpi)\n")
## Figure saved: volcano_plot_LeafxP.png (7 x 7 inches @ 150 dpi)
print(volcano_1panel)

Export Results

Generate LaTeX Tables by Predictor × Regulation

Purpose: Create formatted tables of labeled genes for manuscript inclusion, split by predictor and regulation combination.

Approach: Subset to labeled genes, round numeric values, format for publication, and export separate LaTeX tables for each predictor × regulation combination.

Expected outcome: Multiple LaTeX table files organized by predictor and regulation direction.

# Prepare table data
table_data <- curated %>%
  #filter(!is.na(locus_label) & locus_label != "") %>%
  group_by(predictor, regulation) %>%
  arrange(predictor, regulation, -neglogP) %>%
  mutate(
    logFC = round(logFC, digits = 2),
    neglogP = format(neglogP, digits = 2)
  ) %>%
  select(predictor, regulation, gene, locus_label, desc_merged, neglogP, logFC) %>%
  ungroup()

# Get unique predictor × regulation combinations
pred_reg_combinations <- table_data %>%
  distinct(predictor, regulation) %>%
  arrange(predictor, regulation)

# Generate separate LaTeX table for each combination
for (i in 1:nrow(pred_reg_combinations)) {
  pred <- pred_reg_combinations$predictor[i]
  reg <- pred_reg_combinations$regulation[i]
  
  # Filter data for this combination
  subset_data <- table_data %>%
    filter(predictor == pred, regulation == reg)
  
  # Create filename with sanitized predictor/regulation names
  pred_clean <- gsub(":", "_", pred)
  pred_clean <- gsub("-", "minus", pred_clean)
  reg_clean <- gsub(":", "_", reg)
  reg_clean <- gsub("-", "minus", reg_clean)

  filename <- file.path(paths$tables, paste0("DEG_table_", pred_clean, "_", reg_clean, ".tex"))
  
  # Generate LaTeX table
  subset_data %>%
    select(-predictor, -regulation) %>%
    kbl(
      caption = paste0(
        "Genes with ", tolower(reg), " expression in ", pred, 
        " condition, selected by Mahalanobis distance"
      ),
      format = "latex",
      col.names = c("Gene", "Label", "Description", "-log10(FDR)", "logFC"),
      align = "lllrr",
      booktabs = TRUE
    ) %>%
    kable_classic(full_width = FALSE, latex_options = "HOLD_position") %>%
    save_kable(file = filename, format = "latex")
  
  cat("Saved:", filename, "(", nrow(subset_data), "genes )\n")
}
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_Leaf_Downregulated.tex ( 18 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_Leaf_Upregulated.tex ( 19 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_minusP_Downregulated.tex ( 17 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_minusP_Upregulated.tex ( 17 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_Leaf_minusP_Downregulated.tex ( 13 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_Leaf_minusP_Upregulated.tex ( 11 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_NA_Downregulated.tex ( 0 genes )
## Saved: /Users/fvrodriguez/Library/CloudStorage/GoogleDrive-frodrig4@ncsu.edu/My Drive/repos/inv4m/results/phosphorus_paper/tables/DEG_table_NA_Upregulated.tex ( 0 genes )
# Also generate combined table with all genes
table_data %>%
  kbl(
    caption = paste(
      "All characterized genes selected by Mahalanobis distance outlier detection",
      "across predictors and regulation directions"
    ),
    format = "latex",
    col.names = c("Predictor", "Regulation", "Gene", "Label", "Description",
                  "-log10(FDR)", "logFC"),
    align = "llllrr",
    booktabs = TRUE
  ) %>%
  kable_classic(full_width = FALSE, latex_options = c("HOLD_position", "scale_down")) %>%
  collapse_rows(columns = 1:2, latex_hline = "major") %>%
  save_kable(file = file.path(paths$tables, "DEG_table_all_combined.tex"), format = "latex")

cat("\nSaved: DEG_table_all_combined.tex (", nrow(table_data), "genes total )\n")
## 
## Saved: DEG_table_all_combined.tex ( 122 genes total )

Display Table Previews by Predictor × Regulation

Purpose: Display HTML previews of tables for each predictor × regulation combination in the notebook.

# Display separate tables for each combination
for (i in 1:nrow(pred_reg_combinations)) {
  pred <- pred_reg_combinations$predictor[i]
  reg <- pred_reg_combinations$regulation[i]
  
  # Filter data for this combination
  subset_data <- table_data %>%
    filter(predictor == pred, regulation == reg) %>%
    select(-predictor, -regulation)
  
  # Display header
  cat("\n\n### ", pred, " / ", reg, "\n\n")
  
  # Display table
  tbl <- subset_data %>%
    kbl(
      caption = paste0(nrow(subset_data), " genes"),
      col.names = c("Gene", "Label", "Description", "-log10(FDR)", "logFC")
    ) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "condensed"),
      full_width = FALSE
    )
  
  print(tbl)
}

1 / Downregulated

18 genes
Gene Label Description -log10(FDR) logFC
Zm00001eb152840 pcf7 Transcription factor PCF7 7.5 -1.48
Zm00001eb151160 NA NTF2 domain-containing protein 7.5 -1.15
Zm00001eb076680 sgrl1 Protein STAY-GREEN LIKE, chloroplastic 7.5 -0.95
Zm00001eb038410 ucp4 Mitochondrial uncoupling protein 4 7.5 -0.70
Zm00001eb329970 NA Tyrosine-specific transport protein 7.5 -0.63
Zm00001eb182020 NA protein MAINTENANCE OF PSII UNDER HIGH LIGHT 1 7.5 -0.61
Zm00001eb176730 NA photosynthetic NDH subunit of subcomplex B 1, chloroplastic 7.5 -0.52
Zm00001eb391900 NA Short-chain dehydrogenase TIC 32, chloroplastic 7.4 -0.60
Zm00001eb057540 zmm4 Zea mays MADS4 7.1 -3.40
Zm00001eb154820 NA Choline kinase 7.1 -0.53
Zm00001eb016200 bhlh1 BHLH transcription factor 6.0 -3.62
Zm00001eb364940 plt29 Lipid-transfer protein DIR1 5.2 -2.80
Zm00001eb214750 zmm15 Zea mays MADS-box 15 5.1 -5.04
Zm00001eb320160 alkt1 Alkyl transferase 4.9 -3.82
Zm00001eb169010 ccp18 cysteine protease18 4.0 -2.76
Zm00001eb090330 aatr1 amino acid transporter1 3.8 -3.01
Zm00001eb421180 fp3 Farnesylated protein 3 3.8 -3.23
Zm00001eb411680 glu2 beta-glucosidase2 2.5 -5.12

1 / Upregulated

19 genes
Gene Label Description -log10(FDR) logFC
Zm00001eb297390 hir3 hypersensitive induced reaction3 7.4 0.80
Zm00001eb041700 NA Glycosyltransferase 7.3 1.09
Zm00001eb305330 cyp6 cytochrome P450 7.3 0.90
Zm00001eb037440 bhlh145 bHLH-transcription factor 145 7.0 0.89
Zm00001eb293310 NA DNAJ heat shock N-terminal domain-containing protein 6.7 0.64
Zm00001eb407630 salt1 SalT homolog1 6.5 2.54
Zm00001eb275060 NA NA 6.1 0.73
Zm00001eb098650 trpp2 trehalose-6-phosphate phosphatase2 6.0 1.47
Zm00001eb370960 wrky111 WRKY-transcription factor 111 6.0 0.60
Zm00001eb163980 NA Surfactant protein B containing protein 6.0 0.53
Zm00001eb261620 NA Indole-2-monooxygenase 4.0 2.04
Zm00001eb422900 NA NA 2.8 1.91
Zm00001eb104340 NA MUTL protein homolog 3 2.3 1.96
Zm00001eb169810 NA sphinganine C4-monooxygenase 1 2.2 1.79
Zm00001eb294140 NA NA 2.1 1.90
Zm00001eb002760 NA Cytochrome P450 family 78 subfamily A polypeptide 8 1.8 2.45
Zm00001eb137930 NA 2’-deoxymugineic-acid 2’-dioxygenase 1.6 1.89
Zm00001eb403420 NA ABC-type Co2+ transport system, permease component 1.6 1.81
Zm00001eb054710 NA Chemocyanin 1.5 1.89

2 / Downregulated

17 genes
Gene Label Description -log10(FDR) logFC
Zm00001eb433900 alla1 allantoinase1 6.4 -1.93
Zm00001eb211170 NA Translocase of chloroplast, chloroplastic 5.9 -1.61
Zm00001eb214780 ccp19 cysteine protease19 5.9 -1.95
Zm00001eb070520 bhlh148 bHLH-transcription factor 148 5.8 -2.12
Zm00001eb243180 NA Serine decarboxylase 5.8 -1.74
Zm00001eb377880 NA NA 5.3 -1.63
Zm00001eb114780 cfm3 CRM family member3 4.9 -1.54
Zm00001eb405630 NA C3H transcription factor (Fragment) 4.9 -1.57
Zm00001eb377890 NA SWI/SNF complex component SNF12-like protein 4.8 -1.64
Zm00001eb248820 NA NA 4.7 -1.84
Zm00001eb294690 peamt2 phosphoethanolamine N-methyltransferase 2 3.8 -7.17
Zm00001eb017120 tps8 terpene synthase8 3.3 -4.87
Zm00001eb066620 tut7 Terminal uridylyltransferase 7 2.7 -4.38
Zm00001eb279680 aaap48 amino acid/auxin permease48 2.3 -4.39
Zm00001eb324550 nactf132 NAC-transcription factor 132 2.2 -4.33
Zm00001eb292550 sec14 SEC14 cytosolic factor family protein / phosphoglyceride transfer family protein 1.9 -6.43
Zm00001eb410750 NA 1.4 -4.18

2 / Upregulated

17 genes
Gene Label Description -log10(FDR) logFC
Zm00001eb003820 pilncr1 pi-deficiency-induced long non-coding RNA1 9.0 7.70
Zm00001eb148590 ips1 induced by phosphate starvation1 9.0 7.10
Zm00001eb241920 gpx1 glycerophosphodiester phosphodiesterase1 9.0 6.84
Zm00001eb064450 pap2 purple acid phosphatase2 9.0 4.64
Zm00001eb154650 NA Inorganic pyrophosphatase 1 9.0 3.06
Zm00001eb280120 pfk1 phosphofructose kinase1 9.0 2.58
Zm00001eb063230 plc6 phospholipase C6 9.0 1.90
Zm00001eb313760 NA FLZ-type domain-containing protein 8.9 3.03
Zm00001eb370610 rfk1 Riboflavin kinase 8.9 3.98
Zm00001eb007180 NA Mannose-1-phosphate guanyltransferase alpha 8.8 2.29
Zm00001eb010130 pap19 purple acid phosphatase19 8.8 6.09
Zm00001eb099420 gmps1 GMP synthase 4.3 9.92
Zm00001eb019570 spx7 SPX domain-containing membrane protein7 4.1 8.04
Zm00001eb425050 mdr1 putative multidrug resistance protein 3.6 8.23
Zm00001eb108800 uam1 UDP-arabinopyranose mutase 3.1 8.72
Zm00001eb034810 mgd2 Monogalactosyldiacylglycerol synthase 2.9 11.12
Zm00001eb388800 ltsr1 Low temperature and salt responsive protein 2.3 9.54

3 / Downregulated

13 genes
Gene Label Description -log10(FDR) logFC
Zm00001eb359280 NA NA 5.6 -0.56
Zm00001eb207130 NA Chlorophyll a-b binding protein, chloroplastic 5.4 -1.35
Zm00001eb389720 NA D-fructose-1,6-bisphosphate 1-phosphohydrolase 5.3 -0.81
Zm00001eb070520 bhlh148 bHLH-transcription factor 148 5.1 -0.96
Zm00001eb212520 psad1 photosystem I subunit d1 4.6 -0.62
Zm00001eb179680 NA Chlorophyll a-b binding protein, chloroplastic 4.6 -0.55
Zm00001eb111630 NA Mediator of RNA polymerase II transcription subunit 33A 4.4 -0.60
Zm00001eb362560 ndho1 NADH-plastoquinone oxidoreductase1 4.4 -0.58
Zm00001eb214780 ccp19 cysteine protease19 4.2 -0.82
Zm00001eb071770 mex1 maltose excess protein1 4.0 -0.59
Zm00001eb256120 NA NA 3.8 -1.41
Zm00001eb235450 NA TATA-binding protein-associated factor 2N 3.6 -2.07
Zm00001eb138960 NA NA 2.1 -2.11

3 / Upregulated

11 genes
Gene Label Description -log10(FDR) logFC
Zm00001eb157810 NA Pyruvate kinase 5.6 1.18
Zm00001eb376160 mrpa3 multidrug resistance-associated protein3 5.4 0.64
Zm00001eb063230 plc6 phospholipase C6 4.5 0.56
Zm00001eb144680 NA Ribonuclease T(2) 4.4 0.61
Zm00001eb339870 pld16 phospholipase D16 4.3 0.56
Zm00001eb393060 NA PI-PLC X domain-containing protein 4.0 1.15
Zm00001eb148030 gmp1 mannose-1-phosphate guanylyltransferase1 3.9 0.69
Zm00001eb009430 NA Heptahelical transmembrane protein 4 3.9 0.63
Zm00001eb011050 NA Beta-galactosidase 3.9 0.53
Zm00001eb289800 pah1 phosphatidate phosphatase 1 3.9 0.58
Zm00001eb263160 NA Zinc finger (C3HC4-type RING finger) family protein 2.9 2.16

NA / Downregulated

0 genes
Gene Label Description -log10(FDR) logFC

NA / Upregulated

0 genes
Gene Label Description -log10(FDR) logFC 
cat("\n\nTotal genes across all tables:", nrow(table_data), "\n")

Total genes across all tables: 122

Session Information

sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.6.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US/en_US/en_US/C/en_US/en_US
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] readr_2.1.6      kableExtra_1.4.0 ggtext_0.1.2     ggpubr_0.6.2    
## [5] ggrepel_0.9.6    ggplot2_4.0.1    dplyr_1.1.4      here_1.0.2      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.54          bslib_0.9.0        rstatix_0.7.3     
##  [5] tzdb_0.5.0         vctrs_0.6.5        tools_4.5.1        generics_0.1.4    
##  [9] tibble_3.3.0       pkgconfig_2.0.3    RColorBrewer_1.1-3 S7_0.2.1          
## [13] lifecycle_1.0.4    compiler_4.5.1     farver_2.1.2       stringr_1.6.0     
## [17] textshaping_1.0.4  carData_3.0-5      litedown_0.8       htmltools_0.5.9   
## [21] sass_0.4.10        yaml_2.3.12        Formula_1.2-5      pillar_1.11.1     
## [25] car_3.1-3          jquerylib_0.1.4    tidyr_1.3.1        cachem_1.1.0      
## [29] abind_1.4-8        commonmark_2.0.0   tidyselect_1.2.1   digest_0.6.39     
## [33] stringi_1.8.7      purrr_1.2.0        labeling_0.4.3     cowplot_1.2.0     
## [37] rprojroot_2.1.1    fastmap_1.2.0      grid_4.5.1         cli_3.6.5         
## [41] magrittr_2.0.4     broom_1.0.11       withr_3.0.2        scales_1.4.0      
## [45] backports_1.5.0    rmarkdown_2.30     ggsignif_0.6.4     ragg_1.5.0        
## [49] hms_1.1.4          evaluate_1.0.5     knitr_1.50         viridisLite_0.4.2 
## [53] markdown_2.0       rlang_1.1.6        gridtext_0.1.5     Rcpp_1.1.0        
## [57] glue_1.8.0         xml2_1.5.1         svglite_2.2.2      rstudioapi_0.17.1 
## [61] jsonlite_2.0.0     R6_2.6.1           systemfonts_1.3.1